#!/usr/bin/env python
# coding: utf-8

import sys
import time
import bsddb
import optparse
from collections import defaultdict
from ctypes import cdll, c_double
from os.path import join, abspath, dirname

THPR = cdll.LoadLibrary(abspath(join(dirname(__file__), 'intersect_genes_simple.so')))
THPR.theor_probability_gt.restype    = c_double
THPR.theor_probability_equal.restype = c_double

def calc_thpr(cache, l, n1, n2, n12):
    max_k = min(n1, n2)
    if (n12 <= max_k/2) and (n12 < 10):
        result = 1.0
        for i in xrange(n12):
            x = cache.get(i)
            if x == None:
                x = THPR.theor_probability_equal(l, n1, n2, i)
                cache[i] = x
            result -= x
    else:
        result = 0;
        for i in xrange(n12, max_k+1):
            x = cache.get(i)
            if x == None:
                x = THPR.theor_probability_equal(l, n1, n2, i)
                cache[i] = x
            result += x
    return result
    
def cached_thpr(cache_gt, cache_eq, l, n1, n2, n12):
    probability = cache_gt.get(n12, {}).get(n2)
    if probability == None:
        #probability = THPR.theor_probability_gt(l, n1, n2, n12)
        probability = calc_thpr(cache_eq[n2], l, n1, n2, n12)
        cache_gt[n12][n2] = probability
    return probability


def main(options):
    genomes_num = options.genomes_num
    first_line = open(options.vector_path).readline().strip()
    # check for empty vector
    if first_line:
        new_gene_id, genome_ids_str = first_line.split()
    else:
        return
    vector = genome_ids_str.split(',')

    db = bsddb.btopen(options.matrix_db_path, 'r')
    founded = defaultdict(int)
    for genome_id in vector:
        for gene_id in db[genome_id].split(','):
            founded[gene_id] += 1

    cache_gt = defaultdict(dict)
    cache_eq = defaultdict(dict)    
    n1 = len(vector)
    for gene_id_str, n12 in founded.iteritems():
        gene_id, ref_gene_vector_len = gene_id_str.split(':')
        n2 = int(ref_gene_vector_len)
        if (n12 > float(n1*n2)/genomes_num): # больше среднего
            probability = cached_thpr(cache_gt, cache_eq, genomes_num, n1, n2, n12)
            if probability < 0.01:
                print "%s\t%s\t%g\t%i,%i,%i" % (new_gene_id, gene_id, probability, n1, n2, n12)


if __name__ == "__main__":
    import gc
    gc.disable()

    parser = optparse.OptionParser(
        usage = 'cat results | %prog --gene-map GENES_MAP_FILE --genom-map GENOMES_MAP_FILE',
    )
    parser.add_option(
        '-m',
        dest = 'matrix_db_path',
    )
    parser.add_option(
        '-v',
        dest = 'vector_path',
    )
    parser.add_option(
        '-g',
        dest = 'genomes_num',
        type = 'int',
    )
    opts, args = parser.parse_args()
    if not opts.genomes_num:
        print >> sys.stderr, parser.usage
        sys.exit(1)

    main(opts)

